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Simple practical approach to estimate thermodynamic properties of strongly coupled Yukawa 
systems, in both fluid and solid phases, is presented. The accuracy of the approach is tested 
by extensive comparison with direct computer simulation results (for fluids and solids) and the 
recently proposed shortest-graph method (for solids). Possible applications to other systems of 
softly repulsive particles are briefly discussed. 
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I. INTRODUCTION 

A large class of soft matter systems can be described 
generically as massive charged particles immersed in 
a neutralizing medium of lighter particles Two well- 
known examples are colloidal suspensions^Tii and complex 
(dusty) plasmas.— The interaction between the massive 
particles is a key issue to adequately understand most 
of the phenomena in these systems. Although it can 
be rather complex in real situations, it is often assumed 
that the electrical interactions between the particles can 
be described, in the zero approximation, by the Yukawa 
(also known as the screened Coulomb or Debye-Hiickel) 
repulsive potential 

V (r) = (Q^/r) exp(-r/AD), (1) 

where Q is the particle charge, Ad is the Debye screening 
length associated with the neutralizing medium, and r is 
the distance between a pair of particles. The effect of the 
neutralizing medium on the effective interactions between 
the particles can involve more than only screening. This 
is particularly relevant in complex plasmas, where contin¬ 
uous absorption of plasma electrons and ions on the par¬ 
ticle surface result in inverse-power-law asymptotes of in¬ 
teraction at large interparticle separationsi^^— Plasma 
production and loss processes can also produce a double- 
Yukawa interaction potential characterized by two dif¬ 
ferent screening lengths— Nevertheless, many experi¬ 
mentally observed trends can be already reproduced by 
the simplest model considering point-like particles inter¬ 
acting via the repulsive Yukawa potential m, at least 
qualitatively. For this reason this model, which we will 
refer to as the single component Yukawa system, has re¬ 
ceived a great deal of interest. Another remarkable prop¬ 
erty, which makes the model Yukawa interaction partic¬ 
ularly important for the soft condensed matter research, 
is that by varying Ad an extremely broad range of inter¬ 
action steepness can be explored. The two limiting cases 
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correspond to extremely soft unscreened Coulomb inter¬ 
action (one-component-plasma limit) for Ad —>■ oo and 
extremely steep hard-sphere-like repulsion for Ad —> 0. 

Depending on the strength of the interparticle inter¬ 
action, Yukawa systems can be regarded as weakly or 
strongly coupled. Weak coupling corresponds to the situ¬ 
ation when the energy of the interaction at a characteris¬ 
tic mean interparticle separation A is much smaller than 
the system temperature T (expressed in energy units), 
so that V{A) <^T. In this regime the interaction pro¬ 
vides only small corrections to the thermodynamic quan¬ 
tities of an ideal gas. These corrections can be relatively 
easy evaluated (e.g. by means of the second virial coef¬ 
ficient). In the opposite regime, V{A) I, the inter¬ 
action produce significant (or dominant) contribution to 
the thermodynamic quantities. The particles form either 
the fluid phase or, at even stronger coupling, the solid 
phase. Evaluation of thermodynamic properties becomes 
less trivial than in the limit of weak coupling. The focus 
of the present paper is, therefore, on the strong coupling 
regime. 

Thermodynamic properties and phase diagram of 
strongly coupled Yukawa systems are relatively well in¬ 
vestigated using various computational and analytical 
techniques. Some relevant examples include Monte Carlo 
(MC) and molecular dynamics (MD) numerical simula- 
tions^^— as well as integral equation theoretical stud- 
iesi^Ti^i (see also references therein)—^ In particular, ac¬ 
curate results for the excess internal energy obtained 
in MD and MC simulations have been tabulated in 
Refs. [Ml for a number of state points in the phase dia¬ 
gram of Yukawa systems. However, since the calculation 
of other thermodynamic functions requires knowledge of 
certain integrals and derivatives of the excess energy (see 
examples below), these data are very useful as reference 
points, but less appropriate as a practical tool for es¬ 
timating thermodynamic properties. Semi-empirical fit¬ 
ting formulas2^“— and simplistic approaches^Si^ to esti¬ 
mate thermodynamics of Yukawa systems have been also 
discussed in the literature. These can be useful in cer¬ 
tain situations, but are not very helpful when sufficient 
accuracy is required. 

The purpose of the present paper is to describe simple 
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practical approach to evaluate thermodynamic properties 
of Yukawa systems at strong coupling, both in the fluid 
and solid phases. The approach is based on simple phe¬ 
nomenological arguments, which can be also relevant to a 
wider class of soft repulsive interactions. The accuracy of 
the approach is tested by calculating the compressibility 
factor and comparing this with available as well as new 
results from direct computer simulations. We also use 
the recently proposed shortest-graph methodf^i^ which 
allows to calculate the radial distribution function and, 
hence, pressure in crystals. The relative deviations be¬ 
tween theoretical and simulation results normally do not 
exceed one part in one thousand for dense fluids not too 
far from crystallization and solids not too close to the 
melting point. The largest deviations are documented for 
the solid phase near melting transition and superheated 
solid, but even there they are below 1%. Hence, the pro¬ 
posed approach represents very convenient and accurate 
tool to estimate thermodynamics of Yukawa and other 
related systems in a very broad parameter regime. 


II. THERMODYNAMIC FUNCTIONS OF 
YUKAWA SYSTEMS 


The main thermodynamic quantities of interest here 
are the internal energy U, Helmholtz free energy i^, and 
pressure P of the system of N Yukawa particles, occupy¬ 
ing volume V and having temperature T. These thermo¬ 
dynamic functions are related via^ 



( 2 ) 
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We will use conventional reduced units: u = U/NT, f = 
F/NT, and p = PV/NT. The latter ratio, Z = PV/NT, 
is also known as the compressibility factor. 

The single component Yukawa system is convention¬ 
ally characterized by two dimensionless parameters. The 
first is the coupling parameter, T = Q^/aT, where 
a = (3U/47rY)i/3 is the Wigner-Seitz radius and the 
temperature T is expressed in energy units through¬ 
out this paper. The second is the screening parameter, 
K = o/Ad. It is useful to express reduced thermody¬ 
namic functions in terms of Yukawa system phase state 
variables, k and T. For a fixed number of particles we 
have r oc {aT)~^ oc and k oc a oc . This 

implies 

ar _ r sr _ i r dn _i k 

df~~T' W~~3V' Wr~ ’ W~3V' 


The following equations for the reduced thermodynamic 
functions are then obtained: 
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and 

7 = 1 , 

3ar 3dK 


(5) 


To conclude this section we remind that in this pa¬ 
per we consider exclusively the single component Yukawa 
system, where thermodynamics is completely determined 
by particle-particle interactions and correlations. In real 
systems (e.g. complex plasmas or colloidal suspensions) 
neutralizing medium is normally present to neutralize 
and stabilize like-charged particles. In this case, particle- 
medium interactions also affect thermodynamics. How¬ 
ever, the effect of particle-medium interactions is addi¬ 
tive to that of particle-particles interactions and can be 
easily evaluatedi^ For example, the excess (free) energy 
associated with the presence of neutralizing medium (e.g. 
plasma) is 


^ 3r kF 

/pi - - 2^2 2 • 


( 6 ) 


Note that the contribution of the neutralizing medium 
is negative and dominant at strong coupling, implying 
that the excess energy and pressure of the corresponding 
system are also negative in this regime. For the single 
component Yukawa system these quantities are obviously 
positive. 


III. YUKAWA SOLIDS 

The reduced excess energy of a solid in the harmonic 
approximation is 


u, = M,T+l (7) 

where Mg is the corresponding lattice sum (Madelung 
constant). The reduced Helmholtz free energy of a solid 
is related to the excess energy by the standard integration 
[Eq. dH)], which yields 

/s = Mg(r-r^) + |in(r/F^) + /^, (8) 

where the integration starts from Fm ” the coupling pa¬ 
rameter at the fluid-solid phase transition, and /m = 
/(Fm). According to the Ross melting criterion,— which 
he obtained by reformulating the celebrated Lindemann 
melting law in terms of the statistical-mechanical parti¬ 
tion function 


/^-MgF^+C, (9) 

where C is some constant. In other words, the thermal 
component of the reduced excess free energy remains ap¬ 
proximately constant at melting. Hoover et al^ observed 
that C ~ 6 for several inverse-power-law (IPL) repulsive 
potentials. More recently, Agrawal and Kofke^ docu¬ 
mented some dependence of C on the potential softness 
(C somewhat increases for softer potentials) for a wide 
range of IPL potentials. However, this variation remains 
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FIG. 1: (Color online) Thermal component of the reduced 
excess free energy of Yukawa solids at melting. Red circles 
correspond to the solid bcc lattice, blue crosses denote the 
solid fee phase. Symbols are calculated from the harmonic 
lattice approximation using the data tabulated in Refs. lT^fT^ . 
Solid lines are linear fits to these data points. 


relatively weak, with 5.9 ^ C ^ 6.6 in a very wide range 
of softness investigated. 

We have calculated the values of C at melting of the bcc 
and fee Yukawa solids, in the harmonic lattice approxi¬ 
mation, using the data tabulated in Refs. [TUB In ad¬ 
dition to neglecting anharmonic corrections, we assumed 
that fluid-bcc and fluid-fee transition lines are very close 
to each other in the parameter regime of present inter¬ 
est (for the conhrmation see Fig. 7 of Ref. flit) so that a 
unique value of Fm can be used (this is further justified 
by the fact that Fm appears under the logarithm when 
calculating C). The results are shown in Fig. [TJ Tak¬ 
ing into account some scattering in the data points, it is 
reasonable to apply linear fitting procedure. This yields 


C{k) 


6.55 — 0.13k, (bcc lattice) 
6.57 —0.10k (fee lattice). 


( 10 ) 


The resulting expression for the reduced free energy of 
the solid phase is 

/s(k, f) = M,(K)r + 1 In [r/r„,(K)] + c(k). (ii) 

To use Eq. (HU for practical evaluations we need to spec¬ 
ify the dependence rin(K). A simple expression proposed 
by Vaulina et is suitable for this purpose^ 


rm(K) 


172 exp(aK) 

1 -I- OK -I- ia^K^ ’ 


( 12 ) 


where the constant a = (47r/3)^/^ ~ 1.612 is the ratio 
of the mean interparticle distance A = {V to the 
Wigner-Seitz radius a. 

Equation dS]) complemented by Eqs. (nni-diii) allows 
us to calculate the compressibility of Yukawa solids. We 


TABLE I: The compressibility factor (reduced pressure) for 
the single component Yukawa fee solid. The first two columns 
specify the location of the system state point on the (k, T) 
plane. The third column lists the values of the reduced cou¬ 
pling strength V /Tm, as calculated using Eq. (1121) . (Note that 
the last three rows may correspond to a superheated solid). 
The remaining columns contain the values of Z obtained us¬ 
ing MC simulations by Meijer and Frenkel^^ and tabulated 
in Ref.fiU (Zmc), present theoretical approach (Zpresent), and 
the shortest-graph metho d^^'^^ (.^sg). 


Hi 

F 

F/Fn, 

Zmc 

^present 

^SG 

4.034 

5764.8 

1.44 

41.058 

40.941 

41.197 

4.066 

5719.1 

1.37 

38.927 

38.809 

39.075 

4.096 

5678.3 

1.32 

37.115 

36.992 

37.274 

4.128 

5634.0 

1.26 

35.259 

35.096 

35.399 

4.169 

5577.8 

1.19 

32.969 

32.820 

33.126 

4.208 

5526.8 

1.12 

31.005 

30.870 

31.202 

4.251 

5470.4 

1.05 

28.991 

28.883 

29.196 

4.300 

5408.0 

0.98 

26.919 

26.729 

27.113 

4.344 

5353.2 

0.92 

25.211 

25.000 

25.409 

4.406 

5277.8 

0.84 

22.996 

22.799 

23.268 


compare our theoretical results with the direct MC sim¬ 
ulation of Meijer and FrenkelJ^ In addition, the com¬ 
pressibility has been evaluated using the recently pro¬ 
posed shortest-graph method i^^i^^ In this method the 
pair distribution function g(r) is calculated by assum¬ 
ing that it can be represented as a sum of the Gaussian 
peaks which describe the corresponding partial contribu¬ 
tions of the particles in the crystal. The radial distri¬ 
bution function is then obtained by angular integration, 
g{r) = (l/47r) f dfig(r). The shortest-graph method al¬ 
lows to calculate g(r) for a wide range of interaction 
potentials, crystalline lattices, temperatures, and densi- 
ties^ii^^ The compressibility is obtained from the virial 
equation of state,— 

0 _ AT 

The details of implementation of the shortest-graph 
method are described in Ref. [s^- 

Numerical and theoretical results for the compressibil¬ 
ity of the fee solid phase are summarized in Table ID 
Our present theoretical results are always below the MC 
results, the relative difference increases on approaching 
melting and, in particular, in the superheated regime. 
However, even in this case the difference is below 1%. 
The shortest-graph method demonstrates comparable ac¬ 
curacy, but overestimates the compressibility obtained in 
MC simulations. The documented accuracy is consider¬ 
ably better than (for example) the hard-sphere pertur¬ 
bation theory described in Ref. [H can provide. 

Since the data tabulated previously!^ are limited to a 
relatively narrow range of k around k ~ 4 and are all for 
the fee solid phase, we performed additional comparison 
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TABLE II: The compressibility factor (reduced pressure) for 
the single component Yukawa bcc solid. The first three 
columns are similar to those in Tab||l The remaining columns 
contain the values of Z obtained using our MD simulation 
{Zmd), the shortest-graph metho d^^'^^ (^sg) and the present 
theoretical approach (Zpresent). 


K 

r 

r/Tn, 

Zmd 

.^SG 

■^present 

0.6 

280.6 

1.51 

1090.711 

1090.889 

1090.340 

1.0 

326.1 

1.48 

404.169 

404.297 

403.916 

1.4 

403.2 

1.42 

213.525 

213.637 

213.410 

2.0 

660.1 

1.44 

119.450 

119.447 

119.401 

2.7 

1266.4 

1.41 

73.253 

73.263 

73.290 

4.0 

5755.5 

1.50 

43.258 

43.272 

43.245 


for the regime of smaller k, where Yukawa solid forms a 
stable bcc lattice. The MD simulations have been per¬ 
formed for N = 64000 particles in the canonical (NVT) 
ensemble, using the Langevin thermostat. The numeri¬ 
cal time step is set to where m is the 

particle mass. The cutoff radius for the potential is ISA. 
To obtain the equilibrium state and calculate the com¬ 
pressibility, the simulations have been run for 1.5 x 10® 
steps. 

Comparison between MD simulations, the shortest- 
graph method, and the present theory is summarized 
in Table m The calculations using the shortest-graph 
method yield values extremely close to the “exact” MD 
results. The present analytical approach has almost the 
same accuracy, the agreement with MD results is very im¬ 
pressive. The relative deviations are well below one part 
in one thousand for all phase states. Reasons for these 
deviations can include the approximate character of the 
dependence rin(K) and the neglect of anharmonic correc¬ 
tions to the free energy. These are not easy to estimate. 
Given the simplicity of the present theoretical approach, 
the agreement should be considered as excellent. Thus, it 
can serve as a simple practical tool to evaluate thermody¬ 
namic functions of Yukawa solids when extreme accuracy 
is not required. 


IV. YUKAWA FLUIDS NEAR FREEZING 

The reduced excess energy of the fluid phase can 
be conveniently divided into static and thermal compo- 
nents^i^ii^ 


UR = MflT -I- wth (13) 

The first term corresponds to the static contribution and 
the coefficient Mr is referred to as the fluid Madelung 
constant!^ It can be obtained with the Percus-Yevick 
(PY) radial distribution function of hard spheres in the 
limit 77 = 1 , where 77 is the hard sphere packing frac- 
tionj ^^d® Alternatively it can be evaluated using the ion 
sphere model (ISM), where each particle is placed in the 
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FIG. 2: (Color online) Thermal component of the reduced 
excess energy of Yukawa fluids near freezing versus the re¬ 
duced coupling parameter T/Tm. Symbols correspond to the 
numerical simulations for different values of the screening pa¬ 
rameter re.— The curves correspond to the functional form 
wth = ^(r/Tm)^^® with different values of 5. The value <5 ~ 3.1 
is the most appropriate in the considered near-freezing regime. 


center of the charge neutral Wigner-Seitz spherical cell 
and the energy is then calculated from simple electro¬ 
static consideration.— The result is 


Mr{k) 


k{k + 1 ) 

(re -I- 1) -b (k — ’ 


(14) 


The thermal component of the internal energy exhibits 
quasi-universal scaling for a wide class of soft repulsive 
potentials, as first pointed out by Rosenfeld and Tara- 
zona (RT scaling)— The accuracy of this scaling for 
various model systems has been recently investigated in 
extensive numerical simulations— Previously, a variant 
of the RT scaling has been successively used to obtain 
practical expressions for the internal energy and pres¬ 
sure of Yukawa fluids across coupling regimes— Here 
we modify this approach with the main emphasis on the 
near-freezing regime. 

For strongly coupled Yukawa fluids near freezing the 
RT scaling is illustrated in Fig. [2l where the dependence 
of Uth on F/Fm is plotted for a number of screening pa¬ 
rameters re < 5. Except the strong screening regime with 
re = 5 all other data points have a tendency to group 
around a single quasi-universal curve. Reasonably accu¬ 
rate fits (shown by the curves) can be obtained using the 
functional form 


nth = d(r/F„,)2/5. (15) 

Originally, Rosenfeld suggested^ to use 5 = 3.0, while 
in our previous paper— we used a modified expression, 
nth = 3.2(r/r„,)2/5 — 0.1, in a wider range of coupling. 
Figured] demonstrates that a simple form dTSl) with 5 = 
3.1 is more appropriate near freezing and we use this 
value in calculations. 
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TABLE III: The compressibility factor (reduced pressure) of 
the single component Yukawa fluid (note, however, that the 
first row may correspond to supercooled liquid) at strong cou¬ 
pling. The first three columns are similar to those in Tab. [D 
The remaining columns contain the values of Z obtained us¬ 
ing MC simulations by Meijer and Frenkeli^ and tabulated 
in Ref.fl^ (Zmc), present theoretical approach (Zpresent), and 
our previous approximation (Zprev) from Ref. IdtJ. 


K 

r 

r/r,n 

Zmc 

^present 

^prev 

1.800 

396.9 

1.03 

102.492 

102.498 

102.526 

1.860 

383.9 

0.95 

89.606 

89.542 

89.567 

1.923 

371.4 

0.87 

78.148 

78.123 

78.145 

1.984 

360.0 

0.80 

68.640 

68.618 

68.637 

2.049 

348.6 

0.73 

59.889 

59.880 

59.895 

2.117 

337.5 

0.66 

52.133 

52.138 

52.150 

2.182 

327.3 

0.60 

45.711 

45.699 

45.707 

2.238 

319.2 

0.56 

41.041 

40.998 

41.002 

2.306 

309.7 

0.51 

35.954 

35.904 

35.903 

2.348 

304.2 

0.48 

33.204 

33.189 

33.184 

2.398 

297.9 

0.45 

30.294 

30.258 

30.249 

2.532 

282.1 

0.37 

23.780 

23.760 

23.741 

2.631 

271.5 

0.32 

20.016 

20.017 

19.989 

2.778 

257.1 

0.26 

15.705 

15.724 

15.682 


Integrating the excess energy from F to Fm and sub¬ 
tracting this from the known value of /m(K) we get the 
free energy of Yukawa fluids in the form 


/fl(K,F) =Mfl(K)F- 

y 


Fm(K) [Ms{k) - Mfl(K)] 


Fm(K) 


2/5 


- 1 


-C{k). (16) 


Note that the RT scaling in the form of Eq. (Ha with 
constant S only holds for k < 4 while the fcc-bcc-fluid 
triple point is located near k ~ 4 5 17 , 46,47 lattice 
being thermodynamically favorable at weak screening). 
Therefore, Mhcc{n) and Cbcc(K) should be substituted 
into Eq. (jl 6 l) in the domain of its applicability. 

We calculated the compressibility factors of the 
Yukawa fluid for several strongly coupled state points 
and compared these with the data from Monte Carlo sim¬ 
ulations!^ and with our previous approximation)^ The 
results are summarized in Table Hill Both theoretical ap¬ 
proaches demonstrate excellent agreement with MC sim¬ 
ulations, especially in the vicinity of the fluid-solid phase 
transition. Some advantage of the present construction 
is that the explicit expression for the free energy is avail¬ 
able. The validity of the present approach is merely lim¬ 
ited by the validity of RT scaling, i.e. to the regime of 
strong coupling and weak screening. 
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FIG. 3: (Color online) Approximate static (red solid curve) 
and thermal (blue dashed curve) components of the reduced 
excess energy at freezing of Yukawa fluids. The static contri¬ 
bution dominates for k < 5. 


V. DISCUSSION AND CONCLUSION 

A simple practical approach to estimate thermody¬ 
namic properties of strongly coupled Yukawa systems, 
both in the fluid and solid phases is presented. The ex¬ 
cess free energy is obtained by a standard integration of 
the expressions for the excess internal energy using the 
RT scaling in the dense fluid phase and harmonic lattice 
approximation in the solid phase. The integration starts 
from the freezing/melting point (where the fluid and solid 
free energies are equal by the conventional definition of 
this phase transition in Yukawa systems) and requires 
the knowledge about the location of this point and ex¬ 
cess free energy at this point. This information is avail¬ 
able from earlier numerical simulations. The accuracy of 
the approach is verified by calculating the compressibility 
factors in a wide parameter regime and comparing the re¬ 
sults with those from direct MD and MC computer simu¬ 
lations. The relative deviations are documented to be less 
than 1 % for solids in the near-melting and superheated 
regimes and more than one order of magnitude smaller 
in other regimes investigated. Given the simplicity and 
accuracy of the approach, it represents very convenient 
practical tools to estimate thermodynamic properties of 
Yukawa systems in a very broad range of parameters. 

The important reason behind the success of the present 
approximation is that the static component of the inter¬ 
nal energy dominates over the thermal component for soft 
repulsive interactions. This ensures that some deviations 
from the quasi-universal behavior of Uth for Yukawa fluids 
lead only to very minor corrections of the total internal 
energy. This also guarantees that anharmonic corrections 
to itth in the solid phase do not play very important role. 
To get an idea how soft the repulsion should be, we plot 
in Fig. [3] the static component of the internal energy at 
freezing, MfFm, along with the rough estimate of the 
thermal energy at freezing, itth — 3.0, as functions of the 
screening parameter k. The two curves intersect around 
K ~ 5.5. Thus, the absolute upper boundary of the ap¬ 
proach applicability can be estimated as Atmax — 5. For 













6 


larger k the thermal energy plays significant or domi¬ 
nant role, the RT scaling in its simple form used here is 
no more applicable (see Fig. la, anharmonic corrections 
can play certain role. 

Taking into account certain universal scalings in the 
behavior of simple condensed matter systems^^— we can 
suggest a more general criteria of the applicability of the 
approach described in this paper. We conjecture that 
it should (possibly with some minor modifications) be 
applicable to soft repulsive systems with the generalized 
softness parameter s above certain threshold. The gen¬ 
eralized softness parameter introduced in Ref. [5l| reads 

s = [_l_W'(A)A/W(A)]-\ (17) 

where A characterizes the mean interparticle distance 
and V (r) is the pairwise interaction potential. For exam¬ 
ple, for the inverse-power law (IPL) interactions V(r) oc 
we get the conventional definition s = l/n. Using 
the Yukawa potential and adopting /tmax — 5 we ob¬ 
tain the threshold s* ~ 0.19. This implies for instance, 
that thermodynamics of the IPL systems with n ^ 5 can 
be treated in a manner similar to that applied here to 
Yukawa systems. Indeed, the RT scaling [Eq. (ITSl) with 
constant <5] works well for IPL with n = 4, is slightly less 
appropriate for n = 6, and signihcant deviations are ob¬ 
served for n = 9 and n = 12 (see Fig. 1 of Ref. HI). We 
expect that the present approach can also be helpful in 
connection with the thermodynamics of systems exhibit¬ 
ing water-like anomalies (Gaussian core, exp-6, Hertz and 
related model potentials) in the regime of sufficiently soft 
interaction (high density), but this requires proper ver- 
ihcation. Another interesting question is to which ex¬ 
tent the shortest-graph method is applicable for these 


model systems taking into account their rich solid poly¬ 
morphism. This is also left for future studies. 

The last remark concerns the relevance of the single 
component Yukawa model to real systems, like complex 
(dusty) plasmas or colloidal suspensions. First of all, the 
effect of neutralizing medium has to be properly taken 
into account. This should not constitute a major prob¬ 
lem as we briefly discussed in Sec. HIl Second, the ide¬ 
alized model considered here does not account for some 
important properties of real systems. For example, the 
openness of complex plasma systems is known to be re¬ 
sponsible for some deviations from the model Yukawa 
interaction (as was pointed out in the Introduction) and 
variability of the particle charge. These issues have been 
repeatedly discussed in recent publications and we refer 
the reader to Refs. f^l^ls^ . It is not clear at the moment 
how large the related modihcations to thermodynamic 
properties can be and whether they can be evaluated us¬ 
ing conventional methods. Nevertheless, the present ap¬ 
proach can be an important element in constructing more 
realistic and accurate models for the thermodynamics of 
real systems. 
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